function Step6generate_costShocks_18X31()
% *************************************************************************
% GENERATING COST SHOCKS - 18 X 31: alternative dimension 17M + 1S sectors 31 countries
% *************************************************************************
%       (1) This code collects the following matrices:
%           - PPI
%           - WIOD related matrices: GAMMA, D, ALPHA, M(SIGMA), SECTORWEIGHTS, TOWEIGHTS
%           - Exchange rate
%       (2) Checks country, sector or country-sector alignments across matrices
%       (3) Generates GAMMA matrices for the following counterfactuals:
%           - autarky, - smooth world I, - smooth world II
%           (cost shocks not recovered yet)
%       (4) Generates cost shocks under
%           Cost shock: Baseline with one big services sector: PPI_services = PPI_aggregate, beta_cost = beta_fx = 1, beta_comp = 0
%           Cost shock: Baseline with one big services sector: PPI_services = CPI_services, beta_cost = beta_fx = 1, beta_comp = 0
%       (5) Two sets of aggregated PPI and cost shocks
%           (1) Including Srest
%           (2) Dropping Srest
%       (6) Output is converted to a table & saved in
%           - country-sector level       : 'costShocks_allPeriods_table.xlsx'
%           - aggregated to country level: 'aggregatedCostShocks_allPeriods_table.xlsx'
% *************************************************************************
clear; clc;
tic;

%% Settings
rootfolder      = input('Enter the location for the rootfolder in single quotes (''''): \n');
projectfolder   = [rootfolder '\analysis'];
ppifolder       = [projectfolder '\PPI\MatlabFiles\'];
gammafolder     = [projectfolder '\WIOD\GAMMA\MatlabFiles\'];
alphafolder     = [projectfolder '\WIOD\ALPHA\MatlabFiles\'];
sigmafolder     = [projectfolder '\WIOD\SIGMA\MatlabFiles\'];
sweightsfolder  = [projectfolder '\WIOD\SWEIGHTS\MatlabFiles\'];
toweightsfolder = [projectfolder '\WIOD\TOWEIGHTS\MatlabFiles\'];
fxfolder        = [projectfolder '\FX\MatlabFiles\'];
outputfolder    = [projectfolder '\COSTSHOCKS\18X31\'];
betafolder      = [rootfolder '\data\6_BETA\'];
if exist(outputfolder,'dir') == 0
    mkdir(outputfolder);
end

cd(projectfolder);

% Other input files
beta_file = [betafolder 'sectorspecificbetas.xlsx'];

%% -- (Start looping over the year & month)
fprintf([' * Generating cost shocks (18 X 31) for the year: \n']);

for y = 1995:2011
    yr = num2str(y);
    fprintf([' - ' yr ' \n']);
    
    for m = 1:12
        mo = num2str(m);
        
%         fprintf([' ************************************** \n  * Running year: ' yr ', month: ' mo '  \n ************************************** \n']);
        
        %% Download vectors & matrices
        % one for big services scenarios & full linkages
        % Download Gamma matrix
        GAMMAdata = importdata([gammafolder,'WIOD_gamma_31cty18sec_' yr '.csv']);
        Gamma_text = GAMMAdata.textdata(2:end,1:3);
        GAMMA = GAMMAdata.data;
        
        % List of countries & sectors
        list_country = unique(Gamma_text(:,2));
        n_country = size(list_country,1);
        list_sector = unique(Gamma_text(:,3));
        n_sector = size(list_sector,1);
        n_cs = n_country * n_sector;
        
        % Generate identity & 1 & 0 matrices
        I = eye(size(GAMMA));
        ONE = ones(size(GAMMA));
        ZERO = zeros(size(GAMMA));
        
        % Generate D matrix
        D = diag(1 - sum(GAMMA,2));
        D_inv = D \ I;
        
        % Download PPI vector: srest_agg
        ppi_srest_agg_data = importdata([ppifolder,'PPI_' yr mo '_31cty18sec_srest_agg.csv']);
        ppi_srest_agg_text = ppi_srest_agg_data.textdata(2:end,1:5);
        PPI_sagg = ppi_srest_agg_data.data;
        % -size check
        assert(size(PPI_sagg,1) == n_cs )
        assert(size(PPI_sagg,2) == 1 )
        
        % Download PPI vector: srest_cpi
        ppi_srest_cpi_data = importdata([ppifolder,'PPI_' yr mo '_31cty18sec_srest_cpi.csv']);
        ppi_srest_cpi_text = ppi_srest_cpi_data.textdata(2:end,1:5);
        PPI_scpi = ppi_srest_cpi_data.data;
        % -size check
        assert(size(PPI_scpi,1) == n_cs )
        assert(size(PPI_scpi,2) == 1 )
        
        % Download E matrix - short version
        % E_short is based on country, so not different
        Edata = importdata([fxfolder,'FX_' yr mo '.csv']);
        E_text = Edata.textdata(2:end,1:4);
        E_short = Edata.data;
        % -size check
        assert(size(E_short,1) == n_country )
        assert(size(E_short,2) == 1 )
        
        % Compute E matrix - Kronecker product: E*1's
        E_0 = kron(E_short,ones(n_sector,1));
        
        % Download Alpha matrix
        ALPHAdata = importdata([alphafolder,'WIOD_alpha_31cty18sec_' yr '.csv']);
        Alpha_text = ALPHAdata.textdata(2:end,1:3);
        ALPHA = ALPHAdata.data;
        
        % -row sum check
        try
            assert(round(sum(sum(ALPHA,1))) == size(ALPHA,2));
        catch
            % Exception: SWE sector 19 row sum equals 0 from 2009
            sumZeroAlpha = Alpha_text(sum(ALPHA,1) == 0,1:3);
            exception = [yr {'SWE'} {'19'}];
            assert(sum(strcmp(sumZeroAlpha,exception)) == 3);
            
            ALPHA(507,507) = 1; % put weight 1 to SWE 19 (itself)
            assert(round(sum(sum(ALPHA,1))) == size(ALPHA,2) );
        end
        
        % Generate M (Sigma) matrix
        M = GAMMA ./ (1 - diag(D));
        M (isnan(M)) = 0; % from 2009 we need to fix this
        % -row sum check
        try
            assert(sum(sum(M,2)) == size(M,1));
        catch
            % Exception: SWE sector 19 row sum equals 0
            sumZeroM = Alpha_text(sum(M,2) == 0 |isnan(sum(M,2)) == 1,1:3);
            exception = [yr {'SWE'} {'19'}];
            assert(sum(strcmp(sumZeroM,exception)) == 3);
            assert(round(sum(sum(M,2)),0) == n_cs - 1); % 1 less because one of the row totals = 0
            assert(round(sum(sum(M,2))) == size(M,1) - 1); % 1 less because one of the row totals = 0

            % -Fix -if SWE 19 is NaN
            %{
            sumZeroM = Alpha_text(isnan(sum(M,2)),1:3); 
            sumZeroM = Alpha_text(sum(M,2) == 0),1:3);
            exception = [yr {'SWE'} {'19'}];
            assert(sum(strcmp(sumZeroM,exception)) == 3);
            
            assert(sum(sum(isnan(M),2)) == n_cs);
            M (isnan(M)) = 0;
            assert(sum(sum(M,2)) == size(M,1) - 1); % 1 less because one of the row totals = 0
            %}
        end
        %}
        
        % Download total output weights
        TOWEIGHTSdata = importdata([toweightsfolder,'WIOD_TOweights_31cty18sec_' yr '.csv']);
        TOweights_text = TOWEIGHTSdata.textdata(2:end,1:3);
        TOWEIGHTS = TOWEIGHTSdata.data;
        
%         fprintf('  --- Imported 18 X 31 matrices  \n');
        
        %% Check country - sector alignment across matrices
        
        % Sector list: Gamma versus PPI srest_agg
        gamma_sector = Gamma_text(:,3);
        ppi_sector   = ppi_srest_agg_text(:,5);
        assert(sum(strcmp(ppi_sector,gamma_sector)) == size(gamma_sector,1));
        
        % Sector list: Gamma versus PPI srest_cpi
        gamma_sector = Gamma_text(:,3);
        ppi_sector   = ppi_srest_cpi_text(:,5);
        assert(sum(strcmp(ppi_sector,gamma_sector)) == size(gamma_sector,1));
        
        % Sector list: total output weights versus Gamma
        toweights_sector = TOweights_text(:,3);
        assert(sum(strcmp(toweights_sector,gamma_sector)) == size(gamma_sector,1) );
        
%         fprintf('  --- Sector list: aligned - 18 X 31 \n');
        
        % Country list: Gamma versus PPI srest_agg
        gamma_country = Gamma_text(:,2);
        ppi_country = ppi_srest_agg_text(:,4);
        assert(sum(strcmp(ppi_country,gamma_country)) == size(gamma_country,1));
        
        % Country list: Gamma versus PPI srest_cpi
        gamma_country = Gamma_text(:,2);
        ppi_country = ppi_srest_cpi_text(:,4);
        assert(sum(strcmp(ppi_country,gamma_country)) == size(gamma_country,1));
        
        % Country list: E0 versus Gamma
        E_short_country = E_text(:,4);
        E_country = repelem(E_short_country,n_sector,1);
        assert(sum(strcmp(E_country,gamma_country)) == size(gamma_country,1));
        
        % Country list: total output weights versus Gamma
        toweights_country = TOweights_text(:,2);
        assert(sum(strcmp(toweights_country,gamma_country)) == size(gamma_country,1) );
        
%         fprintf('  --- Country list: aligned - 18 X 31 \n');
        
        %% --- Generate cost shocks ---
        
        %% Cost shock: Baseline with one big services sector: PPI_services = PPI_aggregate, beta_cost = beta_fx = 1, beta_comp = 0
        clear beta_cost beta_fx beta_comp BETA_cost BETA_fx BETA_comp LAMBDA_cost LAMBDA_fx LAMBDA_comp W*
        
        % Betas:
        % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
        beta_cost = 1;
        beta_fx   = 1;
        beta_comp = 0;
        
        % BETA matrices
        % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
        BETA_cost = beta_cost * I;
        BETA_fx   = beta_fx   * I;
        BETA_comp = beta_comp * I;
        
        % Generate LAMBDA
        % (1) '': lambda_cost (2) x: lambda_fx (3) d: lambda_comp
        LAMBDA_cost = I;
        LAMBDA_fx   = I;
        % LAMBDAs (_cost & _fx) generated using the formula are very close to the I -Identity matrix,
        % only when rounded to around the 5th decimal point. Therefore used the Identity matrix directly instead.
        % See below for how close it is.
        
        LAMBDA_comp = ZERO;
        
        % Compute W
        IminusBd_inv  = (I - BETA_comp) \ I ;
        W_first_term  = (LAMBDA_cost + ALPHA' * IminusBd_inv * BETA_comp * BETA_cost * M) \ I;
        W_second_term = (I - BETA_fx - IminusBd_inv * BETA_comp * BETA_fx);
        W_third_term  = (IminusBd_inv * BETA_comp * BETA_fx * M);
        W = W_first_term * (PPI_sagg + ((I - LAMBDA_fx) - ALPHA' * ( W_second_term + W_third_term )) * E_0);
        
        % Check W == PPI
        assert(sum(sum(abs(W - PPI_sagg))) / size (W,1) < 1e-08);
        compare_w_ppi = [W PPI_sagg];
        
        % Compute cost shock, srest: agg
        COST_sagg = D_inv * ((I - IminusBd_inv * BETA_cost * GAMMA) * W + IminusBd_inv * (BETA_fx * (I - D) - BETA_fx * GAMMA ) * E_0);
        simplified_costsrest_agg = D_inv * ((I - GAMMA) * PPI_sagg + (I - D - GAMMA) * E_0);
        
        % Check COST_srest_agg = simplified_cost_srest_agg
        assert(sum(sum(abs(COST_sagg - simplified_costsrest_agg))) / size (COST_sagg,1) < 1e-08);
        compare_costsrestagg_costsrestaggsimplified = [COST_sagg simplified_costsrest_agg];
        
%         fprintf('  --- Generated Cost shocks for baseline with one big services sector: PPI_services = PPI_aggregate - for 18 X 31 \n');
        
        %% Cost shock: Baseline with one big services sector: PPI_services = CPI_services, beta_cost = beta_fx = 1, beta_comp = 0
        clear beta_cost beta_fx beta_comp BETA_cost BETA_fx BETA_comp LAMBDA_cost LAMBDA_fx LAMBDA_comp W*
        
        % Betas:
        % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
        beta_cost = 1;
        beta_fx   = 1;
        beta_comp = 0;
        
        % BETA matrices
        % (1) '': beta_cost (2) x: beta_fx (3) d: beta_comp
        BETA_cost = beta_cost * I;
        BETA_fx   = beta_fx   * I;
        BETA_comp = beta_comp * I;
        
        % Generate LAMBDA
        % (1) '': lambda_cost (2) x: lambda_fx (3) d: lambda_comp
        LAMBDA_cost = I;
        LAMBDA_fx   = I;
        % LAMBDAs (_cost & _fx) generated using the formula are very close to the I -Identity matrix,
        % only when rounded to around the 5th decimal point. Therefore used the Identity matrix directly instead.
        % See below for how close it is.
        
        LAMBDA_comp = ZERO;
        
        % Compute W
        IminusBd_inv  = (I - BETA_comp) \ I ;
        W_first_term  = (LAMBDA_cost + ALPHA' * IminusBd_inv * BETA_comp * BETA_cost * M) \ I;
        W_second_term = (I - BETA_fx - IminusBd_inv * BETA_comp * BETA_fx);
        W_third_term  = (IminusBd_inv * BETA_comp * BETA_fx * M);
        W = W_first_term * (PPI_scpi + ((I - LAMBDA_fx) - ALPHA' * ( W_second_term + W_third_term )) * E_0);
        
        % Check W == PPI
        assert(sum(sum(abs(W - PPI_scpi))) / size (W,1) < 1e-08);
        compare_w_ppi = [W PPI_scpi];
        
        % Compute cost shock, srest: cpi
        COST_scpi = D_inv * ((I - IminusBd_inv * BETA_cost * GAMMA) * W + IminusBd_inv * (BETA_fx * (I - D) - BETA_fx * GAMMA ) * E_0);
        simplified_costsrest_cpi = D_inv * ((I - GAMMA) * PPI_scpi + (I - D - GAMMA) * E_0);
        
        % Check COST_srest_cpi = simplified_cost_srest_cpi
        assert(sum(sum(abs(COST_scpi - simplified_costsrest_cpi))) / size (COST_scpi,1) < 1e-08);
        compare_costsrestcpi_costsrestcpisimplified = [COST_scpi simplified_costsrest_cpi];
        
%         fprintf('  --- Generated Cost shocks for baseline with one big services sector: PPI_services = CPI_services - for 18 X 31 \n');
        
        %% Gather PPI & cost shocks in a table
        headers_textdata = ppi_srest_agg_data.textdata(1,1:5);
        textdata = ppi_srest_agg_data.textdata(2:end,1:5);
        
        % Get the list of all generated cost shocks, and add PPI
        if y == 1995 && m == 1
            numdatanames = who('-regexp','PPI*');
            numdatanames = [numdatanames; who('-regexp','COST*')];
        end
        
        numdata = []; headers_numdata = '';
        for n = 1:length(numdatanames)
            numdata = [numdata eval(numdatanames{n})];
            headers_numdata = [headers_numdata cellstr(numdatanames{n})];
        end
        
        %headers_numdata = ['PPI' 'COST1' ];
        %numdata = [PPI COST1 ];
        
        TABLE = [];
        for i = 1: length(headers_textdata)
            clear T;
            T = table(textdata(:,i),'variableNames',cellstr(headers_textdata{i}));
            TABLE = [TABLE T];
        end
        
        for i = 1: length(headers_numdata)
            clear T;
            T = table(numdata(:,i),'variableNames',cellstr(headers_numdata{i}));
            TABLE = [TABLE T];
        end
        
        %% Generate aggregated PPI & cost shocks at the country level
        % - using 2002 total output weights
        
        % Download total output weights 2002
        TOWEIGHTSdata2002 = importdata([toweightsfolder,'WIOD_TOweights_31cty18sec_2002.csv']);
        TOweights2002_text = TOWEIGHTSdata2002.textdata(2:end,1:3);
        TOWEIGHTS2002 = TOWEIGHTSdata2002.data;
        TOWEIGHTS2002_tomultiply = reshape(TOWEIGHTS2002,[n_sector,n_country]);
        
        % Reshape for multiplying by weighting matrix
        for n = 1:length(numdatanames)
            eval([ numdatanames{n} '_forweighing = transpose(reshape(' numdatanames{n} ',[n_sector,n_country]));']);
        end
        
        % Multiply and take the diagonal: to weigh
        for n = 1:length(numdatanames)
            eval([ numdatanames{n} '_ag = diag(' numdatanames{n} '_forweighing * TOWEIGHTS2002_tomultiply);']);
        end
        
        % --------------------------------
        
        % Drop Srest & compute aggregates
        for n = 1:length(numdatanames)
            eval([ numdatanames{n} '_dropS = ' numdatanames{n} '(strcmp(Gamma_text(:,3),''Srest'') == 0,1);']);
        end
        
        % Download total output weights 2002
        TOWEIGHTSdata2002_17X31 = importdata([toweightsfolder,'WIOD_TOweights_2002.csv']);
        TOweights2002_17X31_text = TOWEIGHTSdata2002_17X31.textdata(2:end,1:3);
        TOWEIGHTS2002_17X31 = TOWEIGHTSdata2002_17X31.data;
        TOWEIGHTS2002_17X31_tomultiply = reshape(TOWEIGHTS2002_17X31,[n_sector-1,n_country]); % 1 sector less: Srest
        
        % Reshape for multiplying by weighting matrix
        for n = 1:length(numdatanames)
            eval([ numdatanames{n} '_dropS_forw = transpose(reshape(' numdatanames{n} '_dropS,[n_sector-1,n_country]));']); % 1 sector less: Srest
        end
        
        % Multiply and take the diagonal: to weigh
        for n = 1:length(numdatanames)
            eval([ numdatanames{n} '_dropS_ag = diag(' numdatanames{n} '_dropS_forw * TOWEIGHTS2002_17X31_tomultiply);']);
        end

        
        %% Gather aggregated PPI & cost shocks in a table
        % Get the list of generated aggregated cost shocks, and add PPI
        numdata_aggregate = []; headers_numdata_aggregate = ''; % initialize
        for n = 1:length(numdatanames)
            numdata_aggregate = [numdata_aggregate eval([numdatanames{n} '_ag']) eval([numdatanames{n} '_dropS_ag'])];
            headers_numdata_aggregate = [headers_numdata_aggregate cellstr([numdatanames{n} '_ag']) cellstr([numdatanames{n} '_dropS_ag'])];
        end
        
        %headers_numdata_aggregate = list2Cell('PPI_aggregate,COST1_aggregate);
        %numdata_aggregate = [PPI_aggregate COST1_aggregate];
        
        % Prepare text data
        selectC = [1:n_sector:n_cs];
        country_text = Gamma_text(selectC,2);
        year_text = cellstr(repmat(yr,size(country_text)));
        month_text = cellstr(repmat(mo,size(country_text)));
        
        % Generate table
        Tcountry = table(country_text,'variableNames',{'country'});
        Tyear    = table(year_text,'variableNames',{'year'});
        Tmonth   = table(month_text,'variableNames',{'month'});
        AGGTABLE = [Tcountry Tyear Tmonth];
        
        for i = 1: length(headers_numdata_aggregate)
            clear T;
            T = table(numdata_aggregate(:,i),'variableNames',cellstr(headers_numdata_aggregate{i}));
            AGGTABLE = [AGGTABLE T];
        end
        
        %% Stack tables for each period in a large table with all periods
        if y == 1995 && m == 1
            TABLE_allP = []; AGGTABLE_allP = []; % initialize
        end
        TABLE_allP = [TABLE_allP; TABLE]; % cs breakdown
        AGGTABLE_allP = [AGGTABLE_allP; AGGTABLE]; % c breakdown
        
        %% -- (End the loop for the month)
    end    
    %% -- (End the loop for the year)
end

%% Export to excel/csv
% Cost shock for each country-sector
output_filename = ['costShocks_allPeriods_table_18X31.xlsx'];
writetable(TABLE_allP,fullfile(outputfolder,output_filename),'WriteRowNames',true);
fprintf(['  --- Cost shocks exported to excel: ' output_filename '  \n']);

% Aggregated cost shock for each country
output_filename = ['aggregatedCostShocks_allPeriods_table_18X31.xlsx']; % The file used in factor analysis
writetable(AGGTABLE_allP,fullfile(outputfolder,output_filename),'WriteRowNames',true);
fprintf(['  --- Aggregated cost shocks exported to excel: ' output_filename '  \n']);

toc
pause(3)
exit
end



